function vector_b = assemble_vector_b_2d(node, elem, function_f, basis_type, test_type)

gauss_weights_tri = [0.068464377671353521259689500766399, 0.10954300427416563401550320122624, 0.068464377671353521259689500766399, 0.061728395061728395061728395061728, 0.098765432098765432098765432098765, 0.061728395061728395061728395061728, 0.0086961161558069665494796751659123, 0.013913785849291147520001565851544, 0.0086961161558069665494796751659123];
gauss_weights_rec = [0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.79012345679012345679012345679012, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864, 0.49382716049382716049382716049383, 0.30864197530864197530864197530864];

vector_b = zeros(size(node,1),1);

psi = reference_basis_points_2d(basis_type, test_type);
ref_trans_tri = reference_basis_points_2d("P1");
ref_trans_rec = reference_basis_points_2d("Q1");

switch basis_type
    case {"P1", "P1b", "P2"}
        for n = 1:size(elem,1)
            E = node(elem(n,:),:);
            V = E(1:3,1:2);
            J_det = ((V(2,1)-V(1,1))*(V(3,2)-V(1,2)) - (V(2,2)-V(1,2))*(V(3,1)-V(1,1)));
            ref_points = [ref_trans_tri{1}(:), ref_trans_tri{2}(:), ref_trans_tri{3}(:)]*V;
            for beta = 1:size(elem,2)
                f = (function_f(ref_points(:,1),ref_points(:,2))' .* psi{beta});
                r = (gauss_weights_tri * f') .* abs(J_det);
                vector_b(elem(n,beta),1) = vector_b(elem(n,beta),1) + r;
            end
        end
    case {"Q1", "Q1b", "Q2"}

    otherwise
        error("Invalid basis type");
end

end